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Previous and present "academic" research aiming at atomic scale understanding is mainly con- 
cerned with the study of individual molecular processes possibly underlying materials science appli- 
cations. In investigations of crystal growth one would for example study the diffusion of adsorbed 
atoms at surfaces, and in the field of heterogeneous catalysis it is the reaction paths of adsorbed 
species that is analyzed. Appealing properties of an individual process are then frequently discussed 
in terms of their direct importance for the envisioned material function, or reciprocally, the function 
of materials is somehow believed to be understandable by essentially one prominent elementary 
process only. What is often overlooked in this approach is that in macroscopic systems of techno- 
logical relevance typically a large number of distinct atomic scale processes take place. Which of 
them are decisive for observable system properties and functions is then not only determined by 
the detailed individual properties of each process alone, but in many, if not most cases also the 
interplay of all processes, i.e. how they act together, plays a crucial role. For a predictive materials 
science modeling with microscopic understanding, a description that treats the statistical interplay 
of a large number of microscopically well-described elementary processes must therefore be applied. 
Modern electronic structure theory methods such as density-functional theory (DFT) have become 
a standard tool for the accurate description of the individual atomic and molecular processes. In 
what follows we discuss the present status of emerging methodologies which attempt to achieve a 
(hopefully seamless) match of DFT with concepts from statistical mechanics or thermodynamics, in 
order to also address the interplay of the various molecular processes. The new quality of, and the 
novel insights that can be gained by, such techniques is illustrated by how they allow the description 
of crystal surfaces in contact with realistic gas-phase environments, which is of critical importance 
for the manufacture and performance of advanced materials such as electronic, magnetic and optical 
devices, sensors, lubricants, catalysts and hard coatings. 



For obtaining an understanding, and for the design, ad- 
vancement or refinement of modern technology that con- 
trols many (most) aspects of our life, a large range of time 
and length scales needs to be described, namely, from 
the electronic (or microscopic/atomistic) to the macro- 
scopic, as illustrated in Figure Obviously, this calls 
for a multi-scale modeling, were corresponding theories 
(i.e. from the electronic, mesoscopic, and macroscopic 
regimes) and their results need to be linked appropriately. 
For each length and time scale regime alone, a number 
of methodologies are well established. It is however, the 
appropriate linking of the methodologies that is only now 
evolving. Conceptually quite challenging in this hierar- 
chy of scales are the transitions from what is often called 
a micro- to a mesoscopic-system description, and from a 
meso- to a macroscopic-system description. Due to the 
rapidly increasing number of particles and possible pro- 
cesses, the former transition is methodologically primar- 
ily characterized by the rapidly increasing importance of 
statistics, while in the latter, the atomic substructure is 
finally discarded in favor of a continuum modeling. In 
this contribution we will concentrate on the micro- to 
mesoscopic-system transition, and correspondingly dis- 
cuss some possibilities of how atomistic electronic struc- 
ture theory can be linked with concepts and techniques 



from statistical mechanics and thermodynamics. 

Since our aim is a materials science modeling that is 
based on understanding, predictive, and applicable to a 
wide range of realistic conditions (e.g. realistic environ- 
mental situations of varying temperatures and pressures), 
this mostly excludes the use of empirical or fitted pa- 
rameters - both at the electronic and at the mesoscopic 
level, as well as in the matching procedure itself. Elec- 
tronic theories that do not rely on such parameters are 
often referred to as first-principles (or in latin: ah initio) 
techniques, and we will maintain this classification also 
for the linked electronic-statistical methods. Correspond- 
ingly, our discussion will mainly (nearly exclusively) fo- 
cus on such ah initio studies, although mentioning some 
other work dealing with important (general) concepts. 
Furthermore, this chapter does not (or only briefly) dis- 
cuss equations; instead the concepts are demonstrated 
(and illustrated) by selected, typical examples. Since 
many (possibly most) aspects of modern material science 
deal with surface or interface phenomena, the examples 
are from this area, addressing in particular surfaces of 
semiconductors, metals, and metal oxides. Apart from 
sketching the present status and achievements, we also 
find it important to mention the difficulties and prob- 
lems (or open challenges) of the discussed approaches. 
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FIG. 1: Schematic presentation of the time and length scales 
relevant for most material science applications. The elemen- 
tary molecular processes, which rule the behavior of a system, 
take place in the so-called "electronic regime" . Their inter- 
play, which frequently determines the functionalities however, 
only develops after meso- and macroscopic lengths or times. 



This can however only be done in a quahtative and rough 

manner, since the problems he mostly in the details, 
the explanations of which are not appropriate for such 
a chapter. 

To understand the elementary processes ruling the ma- 
terials science context, microscopic theories need to ad- 
dress the behavior of electrons and the resulting interac- 
tions between atoms and molecules (often expressed in 
the language of chemical bonds). Electrons move and 
adjust to perturbations on a time scale of femtoseconds 
(1 fs = 10""'^^ s), atoms vibrate on a time scale of picosec- 
onds (1 ps = 10^^^ s), and individual molecular processes 
take place on a length scale of 0.1 nanometer (1 nm = 
10~^m). Because of the central importance of the elec- 
tronic interactions, this time and length scale regime is 
also often called the "electronic regime" , and we will use 
this term here in particular, in order to emphasize the 
subtle difference between ab initio electronic and semi- 
empirical microscopic theories. The former explicitly 
treat the electronic degrees of freedom, while the latter 
already coarse-grain over them and directly describe the 
atomic scale interactions by means of interatomic poten- 
tials. Many materials science applications depend sen- 
sitively on intricate details of bond breaking and mak- 
ing, which on the other hand are often not well (if at 
all) captured by existing semi-empiric classical potential 
schemes. A predictive first-principles modeling as out- 
lined above must therefore be based on a proper descrip- 
tion of molecular processes in the "electronic regime", 
which is much harder to accomplish than just a micro- 
scopic description employing more or less guessed po- 
tentials. In this respect we find it also appropriate to 
distinguish the electronic regime from the currently fre- 
quently cited "nanophysics" (or better "nanometer-scale 



physics"). The latter deals with structures or objects 
of which at least one dimension is in the range 1-100 nm, 
and which due to this confinement exhibit properties that 
are not simply scalable from the ones of larger systems. 
Although the molecular processes in the electronic regime 
operate on a sub-nanometer length scale, they do thus 
not necessarily fall into the "nanophysics" category (un- 
less they take place in a nano-scale system and exhibit 
qualitatively new properties). 

Although already quite involved, the detailed under- 
standing of individual molecular processes arising from 
electronic theories is, however, often still not enough. As 
mentioned above, in many cases the system fmictionali- 
ties are determined by the concerted interplay of many 
elementary processes, not only by the detailed individual 
properties of each process alone. It can for example very 
well be that an individual process exhibits very appeal- 
ing properties for a desired application, yet the process 
may still be irrelevant in practice, because it hardly ever 
occurs within the "full concert" of all possible molecular 
processes. Evaluating this "concert" of elementary pro- 
cesses one obviously has to go beyond separate studies of 
each microscopic process. Taking the interplay into ac- 
count, however, naturally requires the treatment of larger 
system sizes, as well as an averaging over much longer 
time scales. The latter point is especially pronounced, 
since many elementary processes in material sciences are 
activated (i.e. an energy barrier must be overcome) and 
thus rare. This means that the time between consecu- 
tive events can be orders of magnitude longer than the 
actual event time itself. Instead of the above mentioned 
electronic time regime, it can therefore be necessary to 
follow the time evolution of the system up to seconds 
and longer in order to arrive at meaningful conclusions 
concerning the effect of the statistical interplay. Apart 
from the system size, there is thus possibly the need to 
bridge some twelve orders of magnitude in time which 
puts new demands on theories that are to operate in the 
corresponding mesoscopic regime. And also at this level, 
the ab initio approach is much more involved than an em- 
pirical one because it is not possible to simply "lump to- 
gether" several not further specified processes into one ef- 
fective parameter. Each individual elementary step must 
be treated separately, and then combined with all the 
others within an appropriate framework. 

Methodologically, the physics in the electronic regime 
is best described by electronic structure theories, among 
which density-functional theory (Hohenberg and Kohn, 
1964; Kohn and Sham, 1965; Parr and Yang, 1989; Drei- 
zler and Gross, 1990) has become one of the most suc- 
cessful and widespread approaches. Apart from detailed 
information about the electronic structure itself, the typ- 
ical output of such DFT calculations, that is of relevance 
for the present discussion, is the energetics, e.g. total 
energies, as well as the forces acting on the nuclei for 
a given atomic configuration of a microscopically-sized 
system. If this energetic information is provided as func- 
tion of the atomic configuration {R/}, one talks about 
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Information 



Time scale 



Length scale 

< 10^ atoms 

< 10'* atoms 

< 10^ atoms 

< 10^ atoms 
%j 1 /im 

> 10 nm 

> 10 nm 



Density-functional theory 

Ab initio atomistic thermodynamics 

Ab initio molecular dynamics 

Semi-empirical molecular dynamics 

Kinetic Monte Carlo simulations 

Rate equations 

Continuum equations 



Microscopic 
Microscopic 
Microscopic 
Microscopic 
Micro- to mesoscopic 
Averaged 
Macroscopic 



Averaged 
t< 50 ps 
t< I ns 
1 ps Si t ^ 1 hour 
0.1s <t<oo 
Is <t<oo 



TABLE I: The time and length scales typically handled by different theoretical approaches to study chemical reactions and 
crystal growth. 



a potential energy surface (PES) EdRj}). Obviously, a 
(meta)stable atomic configuration corresponds to a (lo- 
cal) minimum of the PES. The forces acting on the given 
atomic configuration are just the local gradient of the 
PES, and the vibrational modes of a (local) minimum 
are given by the local PES curvature around it. Al- 
though DFT mostly does not meet the frequent demand 
for "chemical accuracy" (1 kcal/mol « 0.04eV/atom) in 
the energetics, it is still often sufficiently accurate to al- 
low for the aspired modeling with predictive character. In 
fact, we will see throughout this chapter that error can- 
cellation at the statistical interplay level may give DFT- 
based approaches a much higher accuracy than may be 
expected on the basis of the PES alone. 

With the computed DFT forces it is possible to directly 
follow the motion of the atoms according to Newton's 
laws (Allen and Tildesley, 1997; Frenkel and Smit, 2002). 
With the resulting ab initio molecular dynamics (MD) 
(Car and Parrinello, 1985; Payne et ai, 1992; Galli and 
Pasquarello, 1993; Gross, 1998; Kroes, 1999) only time 
scales up to the order of 50 picoseconds are, however, cur- 
rently accessible. Longer times may e.g. be reached by 
so-called accelerated MD techniques (Voter, Montalenti, 
and Germann, 2002), but for the desired description of 
a truly mesoscopic scale system which treats the statis- 
tical interplay of a large number of elementary processes 
over some seconds or longer, a match or combination of 
DFT with concepts from statistical mechanics or ther- 
modynamics must be found. In the latter approaches, 
bridging of the time scale is achieved by either a suitable 
"coarse-graining" in time (to be specified below) or by 
only considering thcrmodynamically stable states. 

We will discuss how such a description, appropriate 
for a mesoscopic scale system, can be achieved starting 
from electronic structure theory, as well as ensuing con- 
cepts like atomistic thermodynamics, lattice-gas Hamil- 
tonian, equilibrium Monte Carlo simulations, or kinetic 
Monte Carlo simulations. Which of these approaches (or 
a combination) is most suitable depends on the partic- 
ular type of problem. Table U lists the different theo- 
retical approaches and the time and length scales that 
they treat. While the concepts are general, we find it 
instructive to illustrate their power and limitations on 



the basis of a particular issue that is central to the field 
of surface-related studies including applications as im- 
portant as crystal growth and heterogeneous catalysis, 
namely to treat the effect of a finite gas-phase. With 
surfaces forming the interface to the surrounding envi- 
ronment, a critical dependence of their properties on the 
species in this gas-phase, on their partial pressures and 
on the temperature can be intuitively expected (Zangwill, 
1988; Masel, 1996). After all, we recall that for exam- 
ple in our oxygen-rich atmosphere, each atomic site of a 
close-packed crystal surface at room temperature is hit 
by of the order of 10^ O2 molecules per second. That this 
may have profound consequences on the surface structure 
and composition is already highlighted by the every-day 
phenomena of oxide formation, and in humid oxygen-rich 
environments, eventually corrosion with rust and verdi- 
gris as two visible examples (Stampfl et ai, 2002). In 
fact, what is typically called a stable surface structure is 
nothing but the statistical average over all elementary ad- 
sorption processes from, and desorption processes to, the 
surrounding gas-phase. If atoms or molecules of a given 
species adsorb more frequently from the gas-phase than 
they desorb to it, the species' concentration in the surface 
structure will be enriched with time, thus also increas- 
ing the total number of desorption processes. Eventually 
this total number of desorption processes will (averaged 
over time) equal the number of adsorption processes, the 
(average) surface composition and structure will remain 
constant, and the surface has attained its thermodynamic 
equilibrium with the surrounding environment. 

Within this context we may be interested in different 
aspects; for example, on the microscopic level, the first 
goal would be to separately study elementary processes 
such as adsorption and desorption in detail. With DFT 
one could e.g. address the energetics of the binding 
of the gas-phase species to the surface in a variety of 
atomic configurations (Schcffler and Stampfl, 2000), 
and molecular dynamics simulations could shed light on 
the possibly intricate gas-surface dynamics during one 
individual adsorption process (Darling and Holloway, 
1995; Gross, 1998; Kroes, 1999). Already the search for 
the most stable surface structure under given gas-phase 
conditions, however, requires the consideration of the 
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interplay between the elementary processes (of at least 
adsorption and desorption) at the mesoscopic scale. If 
we are only interested in the equilibrated system, i.e. 
when the system has reached its thermodynamic ground 
state, the natural choice would then be to combine DFT 
data with thermodynamic concepts. How this can be 
done will be exemplified in the first part of this chapter. 
On the other hand, the processes altering the surface 
geometry and composition from a known initial state to 
the final ground state can be very slow. And coming back 
to the above example of oxygen-metal interaction, corro- 
sion is a prime example, where such a kinetic hindrance 
significantly slows down (and practically stops) further 
oxidation after a passive oxide film of certain thickness 
has formed at the surface. In such circumstances, a 
thermodynamic description will not be satisfactory and 
one would want to follow the explicit kinetics of the 
surface in the given gas-phase. Then the combination of 
DFT with concepts from statistical mechanics explicitly 
treating the kinetics is required, and we will illustrate 
some corresponding attempts in the last section en- 
titled "First-principles kinetic Monte Carlo simulations" . 

Ab initio atomistic thermodynamics 




FIG. 2: Cartoon sideviews illustrating the effect of an increas- 
ingly oxygen-rich atmosphere on a metal surface. Whereas the 
clean surface prevails in perfect vacuum (left), finite O2 pres- 
sures in the environment also lead to an oxygen-enrichment in 
the solid and its surface. Apart from some bulk dissolved oxy- 
gen, frequently observed stages in this oxidation process com- 
prise (from left to right) on-surface adsorbed O, the forma- 
tion of thin (surface) oxide films, and eventually the complete 
transformation to an ordered bulk oxide compound. Note, 
that all stages can be strongly kinetically-inhibited. It is for 
example, not clear whether the observation of a thin surface 
oxide film means that this is the stable surface composition 
and structure at the given gas-phase pressure and tempera- 
ture, or whether the system has simply not yet attained its 
real equilibrium structure (possibly in form of the full bulk 
oxide). 



Let us at first discuss the matching of electronic struc- 
ture theory data with thermodynamics. Although this 
approach applies "only" to systems in equilibrium, we 
note that at least at not too low temperatures, a sur- 
face is likely to rapidly attain thermodynamic equilib- 
rium with the ambient atmosphere. And even if it has 
not yet equilibrated, at some later stage it will have and 
we can nevertheless learn something by knowing about 
this final state. Thermodynamic considerations also have 
the virtue of requiring comparably less microscopic infor- 
mation, typically only about the minima of the PES and 
the local curvatures around them. As such it is often 
advantageous to first resort to a thermodynamic descrip- 
tion, before embarking upon the more demanding kinetic 
modeling described in the last section. 

The goal of the thermodynamic approach is to use 
the data from electronic structure theory, i.e. the 
information on the PES, to calculate appropriate 
thermodynamic potential functions like the Gibbs free 
energy G (Kaxiras et ai, 1987; Schcfflcr, 1988; Scheffler 
and Dabrowski, 1988; Qian, Martin and Chadi, 1988). 
Once such a quantity is known, one is immediately in 
the position to evaluate macroscopic system properties. 
Of particular relevance for the spatial aspect of our 
multiscale endeavor is further that within a thermody- 
namic description larger systems may readily be divided 
into smaller subsystems that are mutually in equilib- 
rium with each other. Each of the smaller and thus 
potentially simpler subsystems can then first be treated 
separately, and the contact between the subsystems is 
thereafter established by relating their corresponding 
thermodynamic potentials. Such a "divide and conquer" 
type of approach can be especially efficient, if infinite. 



but homogeneous parts of the system like bulk or 
surrounding gas-phase can be separated off (Wang et 
ai, 1998; Wang, Chaka and Scheffler, 2000; Renter and 
Scheffier, 2002; Renter and Scheffler, 2003a, b; Lodzianan 
and N0rskov, 2003). 

Chemical potential plots for surface oxide formation 

How this quite general concept works and what it 
can contribute in practice may be illustrated with the 
case of oxide formation at late transition metal (TM) 
surfaces sketched in Figure |21 (Renter and Scheffler, 
2004a, b). These materials have widespread technologi- 
cal use, for example in the area of oxidation catalysis 
(Ertl, Knozinger and Weitkamp, 1997). Although they 
are likely to form oxidic structures (i.e. ordered oxygen- 
metal compounds) in technologically-relevant high oxy- 
gen pressure environments, it is difficult to address this 
issue at the atomic scale with the corresponding ex- 
perimental techniques of surface science that often re- 
quire Ultra-High Vacuum (UHV) (Woodruff and Delchar, 
1994). Instead of direct, so-called in-situ measurements, 
the surfaces are usually first exposed to a defined oxy- 
gen dosage, and the produced oxygen-enriched surface 
structures are then cooled down and analyzed in UHV. 
Due to the low temperatures, it is hoped that the sur- 
faces do not attain their equilibrium structure in UHV 
during the time of the measurement, and thus provide in- 
formation about the corresponding surface structure at 
higher oxygen pressures. This is, however, not fully cer- 
tain, and it is also not guaranteed that the surface has 
reached its equilibrium structure during the time of oxy- 
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gen exposure at first. Typically, a whole zoo of poten- 
tially kinetically-limited surface structures can be pro- 
duced this way. Even though it can be academically very 
interesting to study all of them in detail, one would still 
like to have some guidance as to which of them would 
ultimately correspond to an equilibrium structure un- 
der which environmental conditions. Furthermore, the 
knowledge of a corresponding, so-called surface phase di- 
agram as a function of in this case the temperature T and 
oxygen pressure po^ can also provide useful information 
to the now surging in-situ techniques, as to which phase 
to expect under which environmental conditions. 

The task for an ab initio atomistic thermodynamic ap- 
proach would therefore be to screen a limited number 
of known (or otherwise provided) oxygen-containing sur- 
face structures, and evaluate which of them turns out 
to be the most stable one under which (T,po2) con- 
ditions (Renter and SchefFler, 2002; Renter and Schef- 
fler, 2003a, b; Lodzianan and N0rskov, 2003). Most sta- 
ble translated into the thermodynamic language meaning 
that the corresponding structure minimizes an appropri- 
ate thermodynamic function, which would in this case 
be the Gibbs free energy of adsorption AG (Li, Stampfl, 
and Scheffler, 2003a,b). In other words, one has to com- 
pute AG as a function of the environmental variables for 
each structural model, and the one with the lowest AG 
is identified as most stable. Obviously this is an indi- 
rect approach, and one of the first limitations of ab initio 
atomistic thermodynamics studies of this kind is thus 
that their reliability is restricted to the structural config- 
urations considered. If the really most stable phase is not 
included in the set of trial structures, the approach will 
not find it, although the obtained phase diagram can well 
give some guidance to what other structures one could 
and should test as well. 

What needs to be computed are all thermodynamic 
potentials entering into the thermodynamic function to 
be minimized. In the present case of the Gibbs free en- 
ergy of adsorption these are for example the Gibbs free 
energies of bulk and surface structural models, as well as 
the chemical potential of the O2 gas-phase. The latter 
may, at the accuracy level necessary for the surface phase 
stability issue, well be approximated by an ideal gas. The 
calculation of the chemical potential fio{T,po2) is then 
elementary and can be found in standard statistical me- 
chanics text books (e.g. McQuarrie, 1976). Required 
input from a microscopic theory like DFT arc properties 
like bond lengths and vibrational frequencies of the gas- 
phase species. Alternatively, the chemical potential may 
be directly obtained from thermochemical tables (StuU 
and Prophet, 1971). Compared to this, the evaluation of 
the Gibbs free energies of the solid bulk and surface is 
more involved. While in principle contributions from to- 
tal energy, vibrational free energy or configurational en- 
tropy have to be calculated (Renter and Scheffler, 2002; 
Reuter and Scheffler, 2003a, b), a key point to notice here 
is that not the absolute Gibbs free energies enter into 
the computation of AG, but only the difference of the 
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FIG. 3: Left panel: Computed Gibbs free energy of adsorp- 
tion AG for the clean Pd(lOO) surface and several oxygen- 
containing surface structures. Depending on the chemical 
potential no of the surrounding gas-phase, cither the clean 
surface or a surface oxide film (labeled here according to its 
twodimensional periodicity as {VEx \/5)R27°), or the infinite 
PdO bulk oxide exhibit the lowest AG and result as the sta- 
ble phase under the corresponding environmental conditions 
(as indicated by the different background shadings). The two 
structures with ordered O adlayers on the surface, p(2 x 2) 
and c(2 x 2), are never most stable, and their frequent ob- 
servation in UHV experiments appears to be the outcome of 
the finite oxygen dosage and kinetic limitations at the low 
temperatures employed. Right panel: The stability range of 
the three phases, evaluated in a) as a function of f^o, plotted 
directly in (T, po2)-spacc. Note the extended stability range 
of the surface oxide compared to the PdO bulk oxide (from 
Reuter and Scheffler, 2004a; Lundgren et al, 2004). 



Gibbs free energies of bulk and surface. This often im- 
plies some error cancellation in the DFT total energies. 
It also leads to quite some (partial) cancellation in the 
free energy contributions like the vibrational energy. In 
a physical picture, it is thus not the effect of the abso- 
lute vibrations that matters for our considerations, but 
only the changes of vibrational modes at the surface as 
compared to the bulk. Under such circumstances it may 
result that the difference between the bulk and surface 
Gibbs free energies is already well approximated by the 
difference of their leading total energy terms, i.e. the 
direct output of the electronic DFT calculations (Reuter 
and Scheffler, 2002). Although this is of course appealing 
from a computational point of view, and one would al- 
ways want to formulate the thermodynamic equations in 
a way that they contain such differences, we stress that it 
is not a general result and needs to be carefully checked 
for every specific system. 

Once the Gibbs free energies of adsorption AG(T, P02 ) 
are calculated for each surface structural model, they 

can be plotted as a function of the environmental con- 
ditions. In fact, under the imposed equilibrium the 
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two-dimensional dependence on T and po^ can be sum- 
marized into a one-dimensional dependence on the gas- 
phase chemical potential noiT^po-i) (Renter and Schef- 
fier, 2002). This is done in Figure for the Pd(lOO) 
surface including, apart from the clean surface, a num- 
ber of previously characterized oxygen-containing surface 
structures. These are two structures with ordered on- 
surface O adsorbate layers of different density [p{2 x 2) 
and c(2 x 2)], a so-called (a/5 x ^/5)R27° surface oxide 
containing one layer of PdO on top of Pd(lOO), and fi- 
nally the infinitely thick PdO bulk oxide (Todorova et 
al, 2003). If we start at very low oxygen chemical poten- 
tial, corresponding to a low oxygen concentration in the 
gas-phase, we expectedly find the clean Pd(lOO) surface 
to yield the lowest AG line, which in fact is used here 
as the reference zero. Upon increasing /io in the gas- 
phase, the Gibbs free energies of adsorption of the other 
oxygen-containing surfaces decrease gradually, however, 
as it becomes more favorable to stabilize such structures 
with more and more oxygen atoms being present in the 
gas-phase. The more oxygen the structural models con- 
tain, the steeper the slope of their AG curves becomes, 
and above a critical /xq we eventually find the surface 
oxide to be more stable than the clean surface. Since the 
PdO bulk oxide contains infinitely many oxygen atoms, 
the slope of its AG line exhibits an "infinite" slope and 
cuts the other lines vertically at A/xq ~ — 0.8eV. For 
any higher oxygen chemical potential in the gas-phase, 
the bulk PdO phase will then always result as most sta- 
ble. 

With the clean surface, the surface and the bulk 
oxide, the thermodynamic analysis yields therefore 
three equilibrium phases for Pd(lGG) depending on the 
chemical potential of the O2 environment. Exploiting 
ideal gas laws, this one dimensional dependence can be 
translated into the physically more intuitive dependence 
on temperature and oxygen pressure. For two fixed 
temperatures, this is also indicated by the resulting 
pressure scales at the top axis of Figure Alterna- 
tively, the stability range of the three phases can be 
directly plotted in (T,po2)"Space, as shown Figure 
Two things are worth noticing. First, in the considered 
thermodynamic equilibrium the two on-surface O adsor- 
bate structures, p{2 x 2) and c(2 x 2), never correspond 
to a stable surface phase, suggesting that their frequent 
observation in UHV experiments is the mere outcome 
of the finite oxygen dosage or kinetic-limitations in the 
afore described preparation procedure. This is indeed an 
apparently unusual result (not yet found for other metal 
surfaces), and, in fact, it is not yet understood. Second, 
the thermodynamic stability range of the recently 
identified surface oxide extends well beyond the one 
of the common PdO bulk oxide, i.e. the surface oxide 
could well be present under environmental conditions 
where the PdO bulk oxide is known to be unstable. Also 
this result is somewhat unexpected, because hitherto 
it had been believed that it is the kinetics (not the 
thermodynamics) that exclusively controls the thickness 



of oxide films at surfaces. The additional stabilization 
of the (\/5 X \/5)R27° surface oxide is attributed to 
the strong coupling of the ultrathin film to the Pd(lOO) 
substrate (Todorova et al., 2003). Similar findings have 
recently been obtained at the Pd(IlI) (Lundgren et ai, 
2002; Reuter and Scheffler, 2004a) and Ag(lll) (Li, 
Stampfl and Scheffler, 2003b; Michaelides et al, 2003) 
surfaces. Interestingly, the low stability of the bulk 
oxide phases of these more noble TMs had hitherto often 
been used as argument against the relevance of oxide 
formation in technological environments like in oxidation 
catalysis (Ertl, Knozinger and Weitkamp, 1997). It 
remains to be seen whether the surface oxide phases 
and their extended stability range, which have recently 
been intensively discussed, will change this common 
perception. 

Chemical potential plots for semiconductor surfaces 

Already in the introduction we had mentioned that the 
concepts discussed here are general and applicable to a 
wide range of problems. To illustrate this, we supplement 
the discussion by an example from the field of semicon- 
ductors, where the concepts of ab initio atomistic ther- 
modynamics had in fact been developed first (Kaxiras et 
al, 1987; Scheffler, 1988; Scheffler and Dabrowski, 1988; 
Qian, Martin and Chadi, 1988). Semiconductor surfaces 
exhibit complex reconstructions, i.e. surface structures 
that differ significantly in their atomic composition and 
geometry from the one that would be obtained by simply 
slice-cutting the bulk crystal (Zangwill, 1988). Knowl- 
edge of the surface atomic structure is, on the other hand, 
a prerequisite to understand and control the surface or 
interface electronic properties, as well as the detailed 
growth characteristics. While the number of possible 
configurations with complex surface unit-cell reconstruc- 
tions is already large, searching for possible structural 
models becomes even more involved for surfaces of com- 
pound semiconductors. In order to minimize the number 
of dangling bonds, the surface may exchange atoms with 
the surrounding gas-phase, which in molecular beam epi- 
taxy (MBE) growth is composed of the substrate species 
at elevated temperatures and varying partial pressures. 
As a consequence of the interaction with this gas-phase, 
the surface stoichiometry may be altered and surface 
atoms be displaced to assume a more favorable bonding 
geometry. The resulting surface structure depends thus 
on the environment, and atomistic thermodynamics may 
again be employed to compare the stability of existing (or 
newly suggested) structural models as a function of the 
conditions in the surrounding gas-phase. The thermo- 
dynamic quantity that is minimized by the most stable 
structure is in this case the surface free energy, which 
in turn depends on the Gibbs free energies of the bulk 
and surface of the compound, as well as on the chemical 
potentials in the gas-phase. The procedure of evaluating 
these quantities goes exactly along the lines described 
above, where in addition, one frequently assumes the sur- 
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FIG. 4: Surface energies for GaAs(OOl) terminations as a 
function of the As chemical potential, /iAs- The thermody- 
namically allowed range of fiAs is bounded by the formation 
of Ga droplets at the surface (As-poor limit at -0.58 eV) and 
the condensation of arsenic at the surface (As-rich limit at 
0.00 eV). The <^(4 x 2) geometry is significantly lower in en- 
ergy than the previously proposed /32(4 x 2) model for the 
c(8 X 2) surface reconstruction observed under As-poor growth 
conditions (from Lee, Moritz and Scheffler, 2000). 



face fringe not only to be in thermodynamic equilibrium 
with the surrounding gas-phase, but also with the under- 
lying compound bulk (Renter and Scheffler, 2002). With 
this additional constraint, the dependence of the surface 
structure and composition on the environment can, even 
for the two component gas-phase in MBE, be discussed 
as a function of the chemical potential of only one of the 
compound species alone. 

Figure ^ shows as an example the dependence on the 
As content in the gas-phase for a number of surface struc- 
tural models of the GaAs(OOl) surface. A reasonable 
lower limit for this content is given, when there is so little 
As2 in the gas-phase that it becomes thermodynamically 
more favorable for the arsenic to leave the compound. 
The resulting GaAs decomposition and formation of Ga 
droplets at the surface denotes the lower limit of As chem- 
ical potentials considered (As-poor limit), while the con- 
densation of arsenic on the surface forms an appropriate 
upper bound (As-rich limit). Depending on the As to 
Ga stoichiometry at the surface, the surface free ener- 
gies of the individual models have either a positive slope 
(As-poor terminations), a negative slope (As-rich termi- 
nations) or remain constant (stoichiometric termination). 
While the detailed atomic geometries behind the consid- 
ered models in Figure ^ are not relevant here, most of 
them may roughly be characterized as different ways of 
forming dimers at the surface in order to reduce the num- 
ber of dangling orbitals (Duke, 1996). In fact, it is this 
general "rule" of dangling bond minimization by dimer 
formation that has hitherto mainly served as inspiration 
in the creation of new structural models for the (001) sur- 
faces of III-V zinc-blende semiconductors, thereby lead- 



ing to some prejudice in the type of structures considered. 
In contrast, the at first theoretically proposed so-called 
C(4 X 2) structure is actuated by the filling of all As dan- 
gling orbitals and emptying of all Ga dangling orbitals, 
as well as a favorable electrostatic (Ewald) interaction 
between the surface atoms (Lee, Moritz and Scheffler, 
2000). The virtue of the atomistic thermodynamic ap- 
proach is now that such a new structural model can be 
directly compared in its stability against all existing ones. 
And indeed, the C(4 x 2) phase was found to be signif- 
icantly more stable than all previously (experimentally) 
proposed reconstructions at low As pressure. 

Returning to the methodological discussion, the results 
shown in Figures 13 and ^nicely summarize the contribu- 
tion that can be made by such ab initio atomistic thermo- 
dynamics considerations. On the other hand, they also 
highlight the limitations. Most prominently, one has to 
be aware that the reliability is restricted to the number 
of considered configurations, or in other words that only 
the stability of exactly those structures plugged in can 
be compared. Had for example the surface oxide struc- 
ture not been known and not explicitly been considered 
in Figure |31 the p{2 x 2) adlayer structure would have 
yielded the lowest Gibbs free energy of adsorption in a 
range of fio intermediate to the stability ranges of the 
clean surface and the bulk oxide, changing the resulting 
surface phase diagram accordingly. Alternatively, it is at 
present not completely clear, whether the (\/5 x \/5)R27° 
structure is really the only surface oxide on Pd(lOO). If 
another yet unknown surface oxide exists and exhibits a 
sufficiently low AG for some oxygen chemical potential, it 
will similarly affect the surface phase diagram, as would 
another novel and hitherto unconsidered surface recon- 
struction with sufficiently low surface free energy in the 
GaAs example. As such, appropriate care should be in 
place when addressing systems where only limited infor- 
mation about surface structures is available. With this 
in mind, even in such systems the atomistic thermody- 
namics approach can still be a particularly valuable tool 
though, since it allows for example to rapidly compare 
the stability of newly devised structural models against 
existing ones. 

In the section entitled "Ab initio lattice-gas Hamil- 
tonian" we will discuss an approach that is able to 
overcome this limitation. This comes unfortunately at 
a significantly higher computational demand, so that it 
has up to now only be used to study simple adsorption 
layers on surfaces. This will then also provide more de- 
tailed insight into the transitions between stable phases. 
In Figures 13 and ^ the transitions are simply drawn 
abrupt, and no reference is made to the finite phase 
coexistence regions that should occur at finite temper- 
atures, i.e. regions in which with changing pressure or 
temperature one phase gradually becomes populated 
and the other one depopulated. That this is not the 
case in the discussed examples has to do with that the 
configurational entropy contribution to the Gibbs free 
energy of the surface phases has been neglected in the 
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corresponding studies (Renter and Scheffler, 2003b). For 
the well-ordered surface structural models considered, 
this contribution is indeed small and will affect only a 
small region close to the phase boundaries. The width 
of this affected phase coexistence region can even be 
estimated, but if more detailed insight into this very 
region is desired, or if disorder becomes more important 
e.g. at more elevated temperatures, then an explicit 
calculation of the configurational entropy contribution 
will become necessary. For this, equilibrium Monte Carlo 
simulations as described below are the method of choice, 
but before we turn to them there is yet another twist to 
chemical potential plots that deserves mentioning. 

"Constrained equilibrium" 

Although a thermodynamic approach can strictly de- 
scribe only the situation where the surface is in equilib- 
rium with the surrounding gas-phase, the idea is that 
it can still give some insight when the system is close to 
thermodynamic equilibrium, or even when it is only close 
to thermodynamic equilibrium with some of the present 
gas-phase species (Renter and SchefSer, 2003a). For such 
situations it can be useful to consider "constrained equi- 
libria" , and one would expect to get some ideas as to 
where in (T,p)-space thermodynamic phases may still 
exist, but also to identify those regions where kinetics 
may control the material function. 

We will discuss heterogeneous catalysis as a prominent 
example. Here, a constant stream of reactants is fed over 
the catalyst surface and the formed products are rapidly 
carried away. If we take the CO oxidation reaction to fur- 
ther specify our example, the surface would be exposed 
to an environment composed of O2 and CO molecules, 
while the produced CO2 desorbs from the catalyst sur- 
face at the technologically employed temperatures and is 
then transported away. Neglecting the presence of the 
CO2, one could therefore model the effect of an O2/CO 
gas-phase on the surface, in order to get some first ideas 
of the structure and composition of the catalyst under 
steady-state operation conditions. Under the assump- 
tion that the adsorption and desorption processes of the 
reactants occur much faster than the CO2 formation re- 
action, the latter would not significantly disturb the aver- 
age surface population, i.e. the surface could be close to 
maintaining its equilibrium with the reactant gas-phase. 
If at all, this equilibrium holds, however, only with each 
gas-phase species separately. Were the latter fully equi- 
librated among each other, too, only the products would 
be present under all environmental conditions of interest. 
It is in fact particularly the high free energy barrier for 
the direct gas-phase reaction that prevents such an equi- 
libration on a reasonable time scale, and necessitates the 
use of a catalyst in the first place. 

The situation that is correspondingly modeled in an 
atomistic thermodynamics approach to heterogeneous 
catalysis is thus a surface in "constrained equilibrium" 
with independent reservoirs representing all reactant gas- 



phase species, namely O2 and CO in the present exam- 
ple (Renter and Scheffler, 2003a). It should immediately 
be stressed though, that such a setup should only be 
viewed as a thought construct to get a first idea about 
the catalyst surface structure in a high-pressure environ- 
ment. Whereas we could write before that the surface 
will sooner or later necessarily equilibrate with the gas- 
phase in the case of a pure O2 atmosphere, this must no 
longer be the case for a "constrained equilibrium" . The 
on-going catalytic reaction at the surface consumes ad- 
sorbed reactant species, i.e. it continuously drives the 
surface populations away from their equilibrium value, 
and even more so in the interesting regions of high cat- 
alytic activity. 

That the "constrained equilibrium" concept can still 
yield valuable insight is nicely exemplified for the CO 
oxidation over a Ru catalyst (Engel and Ertl, 1982). For 
ruthenium the afore described tendency to oxidize un- 
der oxygen-rich environmental conditions is much more 
pronounced than for the above discussed nobler metals 
Pd and Ag (Renter and Scheffler, 2004a). While for 
the latter the relevance of (surface) oxide formation un- 
der the conditions of technological oxidation catalysis is 
still under discussion (Renter and Scheffler, 2004a; Li, 
Stampfl and Schefller, 2003b; Michaehdes et al. 2003; 
Hendriksen, Bobaru and Frenken, 2003), it is by now 
established that a film of bulk-like oxide forms on the 
Ru(OOOl) model catalyst during high-pressure CO oxida- 
tion, and that this RuO2(110) is the active surface for 
the reaction (Over and Muhler, 2003). When evaluating 
its surface structure in "constrained equilibrium" with 
an O2 and CO environment, four different stable phases 
result depending on the gas-phase conditions that are 
now described by the chemical potentials of both reac- 
tants, cf. Figure The phases differ from each other 
in the occupation of two prominent adsorption site types 
exhibited by this surface, called bridge (br) and coordi- 
natively unsaturated (cus) sites. At very low ^j-co, i-e. a 
very low CO concentration in the gas-phase, either only 
the bridge, or bridge and cus sites are occupied by oxy- 
gen depending on the O2 pressure. Under increased CO 
concentration in the gas-phase, both the corresponding 
O^'^/— and the Q^'^ /Q<^^^ phase have to compete with 
CO that would also like to adsorb at the cus sites. And 
eventually the O'^^'/CO'^"'' phase develops. Finally, un- 
der very reducing gas-phase conditions with a lot of CO 
and essentially no oxygen, a completely CO covered sur- 
face results (CO'^'/CO'^"''). Under these conditions the 
RuO2(110) surface can at best be metastable, however, 
as above the white-dotted line in FigureElthe RUO2 bulk 
oxide is already unstable against CO-induced decompo- 
sition. 

With the already described difficulty of operating the 
atomic-resolution techniques of surface science at high 
pressures, the possibility of reliably bridging the so-called 
pressure gap is of key interest in heterogeneous cataly- 
sis research (Ertl, Knozinger, and Weitkamp, 1997; En- 
gel and Ertl, 1982; Ertl, 2002). The hope is that the 
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FIG. 5: Top panel: Top view of the RuO2(110) oxide surface 
explaining the location of the two prominent adsorption sites 
(coordinatively unsaturated, cus, and bridge, br). Also shown 
are side views of the four stable phases present in the phase 
diagram shown below (Ru = light large spheres, O = dark 
medium spheres, C = white small spheres). Bottom panel: 
Surface phase diagram for RuO2(110) in "constrained equi- 
librium" with an oxygen and CO environment. Depending 
on the gas-phase chemical potentials (^o,Mco), br and cus 
sites are either occupied by O or CO, or empty (-), yield- 
ing a total of four different surface phases. For T = 300 K 
and T = 600 K, this dependence is also given in the corre- 
sponding pressure scales. Regions that are expected to be 
particularly strongly affected by phase coexistence or kinetics 
are marked by white hatching (see text). Note that condi- 
tions representative for technological CO oxidation catalysis 
(ambient pressures, 300-600 K) fall exactly into one of these 
ranges (from Renter and Scheffler, 2003a,b). 



atomic-scale understanding gained in experiments with 
some suitably chosen low pressure conditions would also 
be representative of the technological ambient pressure 
situation. Surface phase diagrams like the one shown in 
Figure |S1 could give some valuable guidance in this en- 
deavor. If the {T,po2,Pco) conditions of the low pres- 
sure experiment are chosen such that they lie within the 
stability region of the same surface phase as at high- 



pressures, the same surface structure and composition 
will be present and scalable results may be expected. If, 
however, temperature and pressure are varied in such a 
way, that one crosses from one stability region to another 
one, different surfaces are exposed and there is no rea- 
son to hope for comparable functionality. This would 
e.g. also hold for a naive bridging of the pressure gap by 
simply maintaining a constant partial pressure ratio. 

In fact, the comparability holds not only within the 
regions of the stable phases themselves, but with the 
same argument also for the phase coexistence regions 
along the phase boundaries. The extent of these config- 
urational entropy induced phase coexistence regions has 
been indicated in Figure |S1 by white regions. Although 
as already discussed, the above mentioned approach 
gives no insight into the detailed surface structure under 
these conditions, pronounced fluctuations due to an 
enhanced dynamics of the involved elementary processes 
can generally be expected due to the vicinity of a phase 
transition. Since catalytic activity is based on the same 
dynamics, these regions are therefore likely candidates 
for efhcient catalyst functionality (Renter and Scheffler, 
2003a). And indeed, very high and comparable reaction 
rates have recently been noticed for different environ- 
mental conditions that all lie close to the white region 
between the 0*^70™" and O'^VCO™'' phases. It must 
be stressed, however, that exactly in this region of high 
catalytic activity one would similarly expect the break- 
down of the "constrained equilibrium" assumption of a 
negligible effect of the on-going reaction on the average 
surface structure and stoichiometry. At least everywhere 
in the corresponding hatched regions in Figure |S1 such 
kinetic effects will lead to significant deviations from the 
surface phases obtained within the approach described 
above, even at "infinite" times after steady-state has 
been reached. Atomistic thermodynamics may therefore 
be employed to identify interesting regions in phase 
space. Their surface coverage and structure, i.e. the 
very dynamic behavior, must then however be modeled 
by statistical mechanics explicitly accounting for the 
kinetics, and the corresponding kinetic Monte Carlo sim- 
ulations will be discussed towards the end of the chapter. 

Ab initio lattice- gas Hamiltonian 

The predictive power of the approach discussed in the 
previous sections extends only to the structures that are 
directly considered, i.e., it cannot predict the existence 
of unanticipated geometries or stoichiometrics. It also 
does not explicitly describe coexistence phases or order- 
disorder transitions as configurational entropy is (typ- 
ically) not included. To overcome both of these limita- 
tions, a proper sampling of the whole configuration space 
must be achieved, instead of considering only a set of 
structural models. Modern statistical mechanical meth- 
ods like Monte Carlo (MC) simulations are particularly 
designed to efhciently fulfill this purpose (Frenkel and 
Smit, 2002; Landau and Binder, 2002). The straightfor- 
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FIG. 6: Left: Illustration of some types of lateral interac- 
tions for the case of a two-dimensional adsorbate layer (small 
dark spheres) that can occupy the two distinct threefold hol- 
low sites of a (111) close-packed surface. V^''"^ (n — 1, 2, 3) are 
two-body (or pair) interactions at first, second and third near- 
est neighbor distances of like hollow sites (i.e. fcc-fcc or hcp- 
hcp). V^"° (n = 1, 2, 3) are the three possible three-body (or 
trio) interactions between three atoms in like nearest neighbor 
hollow sites, and V"^^"^'*^''' (^ = 2, 3) represent pair interac- 
tions between atoms that occupy unlike hollow sites (i.e. one 
in fee and the other in hep or vice versa). Right: Exam- 
ple of an adsorbate arrangement from which an expression 
can be obtained for use in solving for interaction parameters. 
The (3 X 3) periodic surface unit-cell is indicated by the large 
darker spheres. The arrows indicate interactions between the 
adatoms. Apart from the obvious first nearest-neighbor inter- 
actions (short arrows), also third nearest-neighbor two-body 
interactions (long arrows) exist, due to the periodic images 
outside of the unit-cell. 



ward matching with electronic structure theories would 
thus be to determine with DFT the energetics of all sys- 
tem configurations generated in the course of the statis- 
tical simulation. Unfortunately, this direct linking is cur- 
rently and also in the foreseeable future computationally 
unfeasible. The exceedingly large configuration spaces of 
most materials science problems require a prohibitively 
large number of free energy evaluations (which can easily 
go beyond 10^ for moderately complex systems). Fur- 
thermore, also disordered configurations must be evalu- 
ated, which in turn is not easy to achieve with the typical 
periodic boundary conditions of DFT supercell calcula- 
tions (SchefHer and Stampfl, 2000). 

With the direct matching impossible, an efficient alter- 
native is to map the real system somehow onto a simpler, 
typically discretized model system, the Hamiltonian of 
which is sufficiently fast to evaluate. This then enables 
us to evaluate the extensive number of free energies re- 
quired by the statistical mechanics. Obvious uncertain- 
ties of this approach are how appropriate the model sys- 
tem represents the real system, and how its parameters 
can be determined from the first-principles calculations. 
The advantage, on the other hand, is that such a de- 
tour via an appropriate ( "coarse-grained" ) model system 
often provides deeper insight and understanding of the 
ruling mechanisms. If the considered problem can be 
described by a lattice defining the possible sites for the 
species in the system, a prominent example for such a 



mapping approach is given by the concept of a lattice- 
gas Hamiltonian (LGH) [or in other languages, an "Ising- 
type model" (de Fontaine, 1994) or a "cluster-expansion" 
(Sanchez, Ducastelle and Gratias, 1984; Zunger, 1994)]. 
Here, any system state is defined by the occupation of 
the sites in the lattice and the total energy of any con- 
figuration is expanded into a sum of discrete interactions 
between these lattice sites. For a one component system 
with only one site type, the LGH would then for example 
read (with obvious generalizations to multi-component, 
multi-site systems): 

pair 

trio 

E E "^"^"'^ + • ■ • ' (1) 

m=l (ijk),„ 

where the site occupation numbers m — or 1 depend 
on whether site I in the lattice is empty or occupied, and 
F is the free energy of an isolated species at this lat- 
tice site, including static and vibrational contributions. 
V^'^" are the two-body (or pair) interaction energies be- 
tween species at mth nearest neighbor sites and is 
the energy due to three-body (or trio) interactions. For- 
mally, higher and higher order interaction terms (quat- 
tro, quinto,...) would follow in this infinite expansion. 
In practice, the series must obviously (and can) be trun- 
cated after a finite number of terms though. Figure El 
illustrates some of these interactions for the case of a 
two-dimensional adsorbate layer that can occupy the two 
distinct threefold hollow sites of a (111) close-packed sur- 
face. In particular, the pair interactions up to third 
nearest neighbor between like and unlike hollow sites are 
shown, as well as three possible trio interactions between 
adsorbates in like sites. 

It is apparent that such a LGH is very general. The 
Hamiltonian can be equally well evaluated for any lattice 
occupation, be it dense or sparse, periodic or disordered. 
And in all cases it merely comprises performing an alge- 
braic sum over a finite number of terms, i.e. it is com- 
putationally very fast. The disadvantage is, on the other 
hand, that for more complex systems with multiple sites 
and several species, the number of interaction terms in 
the expansion increases rapidly. Precisely which of these 
(far-reaching or multi-body) interaction terms need to be 
considered, i.e. where the sum in eq. (1) may be trun- 
cated, and how the interaction energies in these terms 
may be determined, is the really sensitive part of such a 
lattice-gas Hamiltonian approach that must be carefully 
checked. 

The methodology in itself is not new, and traditionally 
the interatomic interactions have often been assumed to 
be just pairwise additive (i.e. higher order terms beyond 
pair interactions were neglected); the interaction ener- 
gies were then obtained by simply fitting to experimental 
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data (see e.g. Piercy, De'Bell, and Pfniir, 1992; Xiong 
et ai, 1997; Zhdanov and Kasemo, 1998). This pro- 
cedure obviously results in "effective parameters" with 
an unsure microscopic basis, "hiding" or "masking" the 
effect and possible importance of three-body (trio) and 
higher-order interactions. This has the consequence that 
while the Hamiltonian may be able to reproduce certain 
specific experimental data to which the parameters were 
fitted, it is questionable and unlikely that it will be gen- 
eral and transferable to calculations of other properties 
of the system. Indeed, the decisive contribution to the 
observed behavior of adparticles by higher-order, many- 
atom interactions has in the meanwhile been pointed out 
by a number of studies (see e.g. Koh and Ehrhch, 1999; 
Osterlund et ai, 1999; Payne et al, 1999). 

As an alternative to this empirical procedure, the lat- 
eral interactions between the particles in the lattice can 
be deduced from detailed DFT calculations, and it is this 
approach in combination with the statistical mechanics 
methods that is of interest for this chapter. The straight- 
forward way to do this is obviously to directly compute 
these interactions as a difference of calculations, in which 
once the involved species are only separately present at 
the corresponding lattice sites, and once all at the same 
time. For the example of a pair interaction between two 
adsorbates at a surface, this would translate into two 
DFT calculations where only either one of the adsorbates 
sits at its lattice site, and one calculation where both 
are present simultaneously. Unfortunately, this type of 
approach is hard to combine with the periodic bound- 
ary conditions that are typically required to describe the 
electronic structure of solids and surfaces (Scheffler and 
Stampfl, 2000). In order to avoid interactions with the 
periodic images of the considered lattice species, huge 
(actually often prohibitively large) supercells would be 
required. A more efficient and intelligent way of address- 
ing the problem is instead to specifically exploit the in- 
teraction with the periodic images. For this, different 
configurations in various (feasible) supercells are com- 
puted with DFT, and the obtained energies expressed in 
terms of the corresponding interatomic interactions. Fig- 
ure illustrates this for the case of two adsorbed atoms 
in a laterally periodic surface unit-cell. Due to this pe- 
riodicity, each atom has images in the neighboring cells. 
Because of these images, each of the atoms in the unit- 
cell experiences not only the obvious pair interaction at 
the first neighbor distance, but also a pair interaction at 
the third neighbor distance (neglecting higher pairwise 
or multi-body interactions for the moment). The com- 
puted DFT binding energy for this configuration i can 
therefore be written as E§^^^'' ^ 2E + 2Vl"'" + 2VP" . 
Doing this for a set of different configurations thus gen- 
erates a system of linear equations that can be solved for 
the interaction energies either by direct inversion (or by 
fitting techniques, if more configurations than interaction 
parameters were determined). 

The crucial aspect in this procedure is the number and 
type of interactions to include in the LGH expansion. 



and the number and type of configurations that are 
computed to determine them. We note that there is no 
a priori way to know at how many, and what type of, 
interactions to terminate the expansion. While there 
are some attempts to automatize this procedure (van de 
Walle and Ceder, 2002), it is probably fair to say that 
the actual implementation remains to date a delicate 
task. Some guidelines to judge on the convergence 
of the constructed Hamiltonian include its ability to 
predict the energies of a number of DFT-computed 
configurations that were not employed in the fit, or that 
it reproduces the correct lowest-energy configurations at 
T OK (so-called "ground-state line"; Zunger, 1994). 

Equilibrium Monte Carlo simulations 

Once an accurate lattice-gas Hamiltonian has been 
constructed, one has at hand a very fast and flexible 
tool to provide the energies of arbitrary system configura- 
tions. This may in turn be used for Monte Carlo simula- 
tions to obtain a good sampling of the available configura- 
tion space, i.e. to determine the partition function of the 
system. An important aspect of modern MC techniques 
is that this sampling is done very efficiently by concen- 
trating on those parts of the configuration space that 
contribute significantly to the latter. The Metropolis al- 
gorithm (Metropolis et ai, 1953), as a famous example 
of such so-called importance sampling schemes, proceeds 
therefore by generating at random new system configu- 
rations. If the new configuration exhibits a lower energy 
than the previous one, it is automatically "accepted" to a 
gradually built-up sequence of configurations. And even 
if the configuration has a higher energy, it still has an 
appropriately Boltzmann weighted probability to make it 
to the sequence. Otherwise it is "rejected" and the last 
configuration copied anew to the sequence. This way, 
the algorithm preferentially samples low energy configu- 
rations, which contribute most to the partition function. 
The acceptance criteria of the Metropolis and of other 
importance sampling schemes are furthermore designed 
in such a way, that they fulfill detailed balance. This 
means that the forward probability of accepting a new 
configuration j from state i is related to the backward 
probability of accepting configuration i from state j by 
the free energy difference of both configurations. Taking 
averages of system observables over the thus generated 
configuration sequences yields then their correct thermo- 
dynamic average for the considered ensemble. Techni- 
cal issues regard finally how new trial configurations are 
generated, or how long and in what system size the sim- 
ulation must be run in order to obtain good statistical 
averages (Frenkel and Smit, 2002; Landau and Binder, 
2002). 

The kind of insights that can be gained by such a first- 
principles LGH -|- MC approach is nicely exemplified by 
the problem of on-surface adsorption at a close-packed 
surface, when the latter is in equilibrium with a surround- 
ing gas-phase. If this environment consists of oxygen. 
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FIG. 7: Phase diagram for O/Ru(0001) as obtained using 
the ab initio lattice-gas Hamiltonian approach in combination 
with MC calculations. The triangles indicate first order tran- 
sitions and the circles second order transitions. The identified 
ordered structures are labeled as: (2x2)-0 (A), (2xl)-0 (B), 
{VSx VS)R30° (C), (2 X 2)-30 (D), and disordered lattice-gas 
(I.g.). (From McEwen, Payne and Stampfl, 2002). 



this would e.g. contribute to the understanding of one 
of the early oxidation stages sketched in Figure |21 What 
would be of interest is for instance to know how much 
oxygen is adsorbed at the surface given a certain tem- 
perature and pressure in the gas-phase, and whether the 
adsorbate forms ordered or disordered phases. As out- 
lined above, the approach proceeds by first determining 
a LGH from a number of DFT-computed ordered adsor- 
bate configurations. This is followed by grand-canonical 
MC simulations, in which new trial system configurations 
are generated by randomly adding or removing adsor- 
bates from the lattice positions and where the energies 
of these configurations are provided by the LGH. Eval- 
uating appropriate order parameters that check on pre- 
vailing lateral periodicities in the generated sequence of 
configurations, one may finally plot the phase diagram, 
i.e. what phase exists under which (r,p)-conditions (or 
equivalently (T, /i)-conditions) in the gas-phase. 

The result of one of the first studies of this kind is 
shown in Figure[71for the system O/Ru(0001). The em- 
ployed lattice-gas Hamiltonian comprised two types of 
adsorption sites, namely the hep and fee hollows, lateral 
pair interactions up to third neighbor and three types 
of trio interactions between like and unlike sites, thus 
amounting to a total of fifteen independent interaction 
parameters. At low temperature, the simulations yield a 
number of ordered phases corresponding to different pe- 
riodicities and oxygen coverages. Two of these ordered 
phases had already been reported experimentally at the 
time the work was carried out. The prediction of two 
new (higher coverage) periodic structures, namely a 3/4 
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FIG. 8: Theoretical (left panel) and experimental (right 
panel) temperature programmed desorption curves. Each 
curve shows the rate of oxygen molecules that desorb from the 
Ru(OOOl) surface as a function of temperature, when the sys- 
tem is prepared with a given initial oxygen coverage 9 ranging 
from 0.1 to 1 monolayer. The first-principles LGH employed 
in the calculations is exactly the same as the one underlying 
the phase diagram of Figure Q (from Stampfi et al, 1999a,b). 



and a 1 monolayer phase, has in the meanwhile been 
confirmed by various experimental studies. This exam- 
ple thus demonstrates the predictive nature of the first- 
principles approach, and the stimulating and synergetic 
interplay between theory and experiment. It is also worth 
pointing out that these new phases and their coexistence 
in certain coverage regions were not obtained in early 
MC calculations of this system based on an empirical 
LGH, which was determined by simply fitting a minimal 
number of pair interactions to the then available experi- 
mental phase diagram (Piercy, De'Bell and Pfniir, 1992). 
We also like to stress the superior transferability of the 
first-principles interaction parameters. As an example 
we name simulations of temperature programmed des- 
orption (TPD) spectra, which can among other possibil- 
ities be obtained by combining the LGH with a transfer- 
matrix approach and kinetic rate equations (Kreuzer and 
Payne, 2000). Figure |H1 shows the result obtained with 
exactly the same LGH that also underlies the phase dia- 
gram of Figure [7| Although empirical fits of TPD spec- 
tra may give better agreement between calculated and 
experimental results, we note that the agreement visible 
in Figure IHl is in fact quite good. The advantage, on the 
other hand, is that no empirical parameters were used 
in the LGH, which allows to unambiguously trace back 
the TPD features to lateral interactions with well defined 
microscopic meaning. 

The results summarized in Figure [7| serve also quite 
well to illustrate the already mentioned differences be- 
tween the initially described chemical potential plots and 
the LGH -I- MC method. In the first approach the sta- 
bility of a fixed set of configurations is compared in order 
to arrive at the phase diagram. For the O/Ru(0001) sys- 
tem, the likely choice at the time would have been just the 
two experimentally known ordered phases, 0(2 x 2) and 
0(2 X 1). The stability region of the prior phase, bounded 
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at lower chemical potentials by the clean surface and at 
higher chemical potentials by the 0(2 x 1) phase, then 
comes out just as much as in Figure □ This stabihty 
range will be independent of temperature, however, i.e. 
there is no order-disorder transition at higher tempera- 
ture due to the neglect of configurational entropy. More 
importantly, since the two new higher-coverage phases 
would not have been explicitly considered, the stability 
of the 0(2 X 1) phase would falsely extend over the whole 
higher chemical potential range. While this emphasizes 
the deeper insight and increased predictive power that is 
achieved by the proper sampling of configuration space 
in the LGH -|- MC technique, one must also recognize 
that the computational cost of the latter is significantly 
higher. It is in particular straightforward to directly com- 
pare the stability of qualitatively different geometries like 
the on-surface adsorption and the surface oxide phases 
in Figure |21 in a chemical potential plot. Setting up a 
lattice-gas Hamiltonian that would equally describe both 
systems, on the other hand, is far from trival. Even if 
it were feasible to find a generalized lattice that would 
be able to treat all system states, disentangling and de- 
termining the manifold of interaction energies in such a 
lattice will be very involved. The required discretiza- 
tion of the real system, i.e. the mapping onto a lattice, 
is therefore to date the major limitation of the LGH + 
MC technique - be it applied to two-dimensional pure 
surface systems or even worse to three-dimensional prob- 
lems addressing a surface fringe of finite width. Still, it 
is also precisely this mapping and the resulting very fast 
analysis of the properties of the LGH that allows for an 
extensive and reliable sampling of the configuration space 
of complex systems that is hitherto unparalleled by other 
approaches. 

Having highlighted the importance of this sampling for 
the determination of unanticipated new ordered phases 
at lower temperatures, the final example in this sec- 
tion illustrates specifically the decisive role it plays also 
for the simulation and understanding of order-disorder 
transitions at elevated temperatures. A particularly in- 
triguing transition of this kind is observed for Na on 
Al(OOl). The interest in such alkali metal adsorption 
systems has been intense, especially since in the early 
1990's it was found (first for Na on Al(lll) and then on 
Al(lOO)) that the alkali metal atoms may kick-out sur- 
face Al atoms and adsorb substitutionally. This was in 
sharp contrast to the generally accepted understanding 
of the time, which was that alkali-metal atoms adsorb 
in the highest coordinated on-surface hollow site, and 
cause little disturbance to a close-packed metal surface 
(Stampfl and Scheffler, 1995; Adams, 1996; Diehl and Mc 
Grath, 1997). For the specific system Na on Al(OOl) at 
a coverage of 0.2 monolayers, a reversible phase transi- 
tion is observed in low energy electron diffraction experi- 
ments at r = 240 K. Below this temperature, an ordered 
(\/5 X \/5)R27° structure forms, where the Na atoms oc- 
cupy surface substitutional sites. At temperatures above 
240 K on the other hand, the Na atoms, still in the same 
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FIG. 9: Top panel: Calculated logarithm of the density of 
configuration states, g{E), for Na on Al(lOO) at a coverage 
of 0.2 monolayers. Lower panels: Internal and free energy of 
this system as derived from g{E) (from Borg et al, 2004). 



substitutional sites, form a disordered arrangement in the 
surface. 

Using the ab initio LGH -f- MC approach the or- 
dered phase and the disorder transition can be success- 
fully reproduced. Pair interactions up to the ninth near- 
est neighbor and two different trio interactions were in- 
cluded in the LGH expansion. To specifically identify 
the crucial role played by configurational entropy in the 
temperature induced order-disorder transition, a specific 
MC algorithm proposed by Wang and Landau (Wang 
and Landau, 2001) was employed. In contrast to the 
above outlined Metropolis algorithm, this scheme af- 
fords an explicit calculation of the density of configu- 
ration states, g{E), i.e. the number of system configura- 
tions with a certain energy E. This quantity provides in 
turn all major thermodynamic functions, e.g., the canon- 
ical distribution at a given temperature, g{E)e~-^^''^'^ , 
the free energy, F{T) = -kBTlniJ^E 9(^)6-^/''^'^) = 
—ksT ln{Z), where Z is the partition function, the in- 
ternal energy, U{T) = [J2e ^gi^)'^''^^''"'^]/^ ^ and the 
entropy 5"= {U -F)/T. 

Figure shows the calculated density of configuration 
states g{E), together with the internal and free energy 
derived from it. In the latter two quantities, the 
abrupt change corresponding to the first-order phase 
transition obtained at 210 K can nicely be discerned. 
In particular, the free energy decreases notably with 
increasing temperature. The reason for this is clearly the 
entropic contribution (difference in the free and internal 
energies), the magnitude of which suddenly increases 
at the transition temperature and continues to increase 
steadily thereafter. Taking this configurational entropy 
into account is therefore (and obviously) the crucial 
aspect in the simulation and understanding of this 
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order-disorder phase transition, and only the LGH+MC 
approach with its proper samphng of configuration space 
can provide it. What the approach does not yield, on 
the other hand, is how the phase transition actually 
takes place microscopically, i.e. how the substitutional 
Na atoms move their positions by necessarily displacing 
surface Al atoms, and on what time scale (with what 
kinetic hindrance) this all happens. For this, one 
necessarily needs to go beyond a thermodynamic de- 
scription, and explicitly follow the kinetics of the system 
over time, which will be the topic of the following section. 

First-principles kinetic Monte Carlo simula- 
tions 
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Up to now we had discussed how equilibrium Monte 
Carlo simulations can be used to explicitly evaluate the 
partition function, in order to arrive at surface phase dia- 
grams as function of temperature and partial pressures of 
the surrounding gas-phase. For this, statistical averages 
over a sequence of appropriately sampled configurations 
were taken, and it is appealing to also connect some time 
evolution to this sequence of generated configurations 
(MC steps). In fact, quite a number of non-equilibrium 
problems may already be tackled on the basis of this un- 
calibrated "MC time" (Landau and Binder, 2002). The 
reason why this does not work in general is twofold. First, 
equilibrium MC is designed to achieve an optimum sam- 
pling within configurational phase space. As such, also 
MC moves that are unphysical like a particle hop from 
an occupied site to an unoccupied one, hundreds of lat- 
tice spacings away may be allowed, if they help to ob- 
tain an efficient sampling of the relevant configurations. 
The remedy for this obstacle is straightforward, though, 
as one only needs to restrict the possible MC moves to 
"physical" elementary processes. The second reason on 
the other hand is more involved, as it has to do with 
the probabilities with which the individual events are ex- 
ecuted. In equilibrium MC the forward and backward 
acceptance probabilities of time-reversed processes like 
hops back and forth between two sites only have to ful- 
fill the detailed balance criterion, and this is not enough 
to establish a proper relationship between MC time and 
"real time" (Kang and Weinberg, 1995). 

In kinetic Monte Carlo simulations (kMC) a proper 
relationship between MC time and real time is achieved 
by interpreting the Monte Carlo process as providing a 
numerical solution to the Markovian master equation 
describing the dynamic system evolution (Bortz, Kalos, 
and Lebowitz, 1975; Gillespie, 1976; Voter, 1986; Kang 
and Weinberg, 1989; Fichthorn and Weinberg, 1991). 
The simulation itself still looks superficially similar to 
equilibrium Monte Carlo in that a sequence of config- 
urations is generated using random numbers. At each 
configuration, however, all possible elementary processes 
and the rates with which they occur are evaluated. Ap- 
propriately weighted by these different rates one of the 
possible processes is then executed randomly to achieve 



FIG. 10: Flow diagram illustrating the basic steps in a ki- 
netic Monte Carlo simulation, i) Loop over all lattice sites of 
the system and determine the atomic processes that are possi- 
ble for the current system configuration. Calculate or lookup 
the corresponding rates, ii) Generate two random numbers, 
iii) advance the system according to the process selected by 
the first random number (this could e.g. be moving an atom 
from one lattice site to a neighboring one, if the correspond- 
ing difi'usion process was selected), iv) Increment the clock 
according to the rates and the second random number, as pre- 
scribed by an ensemble of Poisson processes, and v) start all 
over or stop, if a sufficiently long time span has been simu- 
lated. 



the new system configuration, as sketched in Figure 
[TUl This way, the kMC algorithm effectively simulates 
stochastic processes described by a Poisson distribution, 
and a direct and unambiguous relationship between 
kMC time and real time can be established (Fichthorn 
and Weinberg, 1991). Not only does this open the 
door to a treatment of the kinetics of non-equilibrium 
problems. It also does so very efficiently, since the time 
evolution is actually coarse-grained to the really decisive 
rare events, passing over the irrelevant short-time 
dynamics. Time scales of the order of seconds or longer 
for mesoscopically-sized systems are therefore readily 
accessible by kMC simulations (Voter, Montalenti, and 
Germann, 2002). 

Insights from MD, MC and kMC 

To further clarify the different insights provided by 
molecular dynamics, equilibrium and kinetic Monte 
Carlo simulations, consider the simple, but typical rare 
event type model system shown in Figure 1111 An iso- 
lated adsorbate vibrates at finite temperature T with 
a frequency on the picosecond time scale and diffuses 
about every microsecond between two neighboring sites 
of different stability. In terms of a PES, this situation is 
described by two stable minima of different depths sepa- 
rated by a sizable barrier. Starting with the particle in 
any of the two sites, a MD simulation would follow the 
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FIG. 11: Schematic potential energy surface (PES) represent- 
ing the thermal diffusion of an isolated adsorbate between two 
stable lattice sites A and B of different stability. A MD sim- 
ulation would explicitly follow the dynamics of the vibrations 
around a minimum, and is thus inefficient to address the rare 
diffusion events happening on a much longer time scale. Equi- 
librium Monte Carlo simulations provide information about 
the average thermal occupation of the two sites, <A^>, based 
on the depth of the two PES minima {Ea and Eb)- Kinetic 
Monte Carlo simulations follow the "coarse-grained" time evo- 
lution of the system, N{t), employing the rates for the diffu- 
sion events between the minima (rA-»B, ^b^a). For this, PES 
information not only about the minima, but also about the 
barrier height at the transition state (TS) between initial and 
final state is required (A_Ba, A_Bb). 



thermal motion of the adsorbate in detail. In order to 
do this accurately, timesteps in the femtosecond range 
are required. Before the first diffusion event can be ob- 
served at all, of the order of 10^ time steps have therefore 
to be calculated first, in which the particle does nothing 
but just vibrate around the stable minimum. Computa- 
tionally this is unfeasible for any but the simplest model 
systems, and even if it were feasible it would obviously 
not be an efficient tool to study the long-term time evo- 
lution of this system. 

For Monte Carlo simulations on the other hand, the 
system first has to be mapped onto a lattice. This is un- 
problematic for the present model and results in two pos- 
sible system states with the particle being in one or the 
other minimum. Equilibrium Monte Carlo provides then 
only time-averaged information about the equilibrated 
system. For this, a sequence of configurations with the 
system in either of the two system states is generated, 
and considering the higher stability of one of the min- 
ima, appropriately more configurations with the system 
in this state are sampled. When taking the average, one 
arrives at the obvious result that the particle is with a 
certain higher (Boltzmann-weighted) probability in the 
lower minimum than in the higher one. 

Real information on the long term time-evolution of 
the system, i.e. focusing on the rare diffusion events, 
is finally provided by kinetic Monte Carlo simulations. 
For this, first the two rates of the diffusion events from 
one system state to the other and vice versa have to be 



known. We will describe below that they can be ob- 
tained from knowledge of the barrier between the two 
states and the vibrational properties of the particle in 
the minima and at the barrier, i.e. from the local cur- 
vatures. A lot more information on the PES is therefore 
required for a kMC simulation than for equilibrium MC, 
which only needs input about the PES minima. Once 
the rates are known, a kMC simulation starting from 
any arbitrary system configuration will first evaluate all 
possible processes and their rates and then execute one 
of them with appropriate probability. In the present ex- 
ample this list of events is trivial, since with the particle 
in either minimum only the diffusion to the other mini- 
mum is possible. When the event is executed, on average 
the time (rate)~^ has passed and the clock is advanced 
accordingly. Note that as described initially, the rare dif- 
fusion events happen on a time scale of microseconds, i.e. 
with only one executed event the system time will be di- 
rectly incremented by this amount. In other words, the 
time is coarse-grained to the rare event time, and all the 
short-time dynamics (corresponding in the present case 
to the picosecond vibrations around the minimum) are 
efficiently contained in the process rate itself. 

Since the barrier seen by the particle when in the shal- 
lower minimum is lower than when in the deeper one, 
cf. Figure 1111 the rate to jump into the deeper mini- 
mum will correspondingly be higher than the one for the 
backwards jump. Generating the sequence of configu- 
rations, each time more time will therefore have passed 
after a diffusion event from deep to shallow compared 
to the reverse process. When taking a long-time average, 
describing then the equilibrated system, one therefore ar- 
rives necessarily at the result that the particle is on aver- 
age longer in the lower minimum than in the higher one. 
This is identical to the result provided by equilibrium 
Monte Carlo, and if only this information is required, 
the latter technique would most often be the much more 
efficient way to obtain it. KMC, on the other hand, has 
the additional advantage of shedding light on the detailed 
time-evolution itself, and can in particular also follow the 
explicit kinetics of systems that are not (or not yet) in 
thermal equilibrium. 

From the discussion of this simple model system, it is 
clear that the key ingredients of a kMC simulation are 
the analysis and identification of all possibly relevant 
elementary processes and the determination of the 
associated rates. Once this is known, the coarse graining 
in time achieved in kMC immediately allows to follow 
the time evolution and the statistical occurrence and 
interplay of the molecular processes of mesoscopically 
sized systems up to seconds or longer. As such it is 
currently the most efficient approach to study long 
time and larger length scales, while still providing 
atomistic information. In its original development, kMC 
was exclusively applied to simplified model systems, 
employing a few processes with guessed or fitted rates 
(see e.g. Kang and Weinberg, 1995). The new aspect 
brought into play by so-called first-principles kMC 
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FIG. 12: Calculated DFT-PES of a CO oxidation reaction 
process at the RuO2(110) model catalyst surface. The high- 
dimensional PES is projected onto two reaction coordinates, 
representing two lateral coordinates of the adsorbed O'^"'' and 
QQcus Figure 0. The energy zero corresponds to the 

initial state at (0.00 A, 3.12 A), and the transition state is at 
the saddle point of the PES, yielding a barrier of 0.89 eV. De- 
tails of the corresponding transition state geometry are shown 
in the inset. Ru — light, large spheres, O — dark, medium 
spheres, and C = small, white spheres (only the atoms ly- 
ing in the reaction plane itself are drawn as threedimensional 
spheres). (From Renter and SchefHer, 2003b). 

simulations (Ruggerone, Ratsch and Scheffler, 1997; 
Ratsch, Ruggerone and Scheffler, 1998) is that these 
rates and the processes are directly provided from 
electronic structure theory calculations, i.e. that the 
parameters fed into the kMC simulation have a clear 
microscopic meaning. 

Getting the processes and their rates 

For the rare event type molecular processes mostly 
encountered in the surface science context, an efficient 
and reliable way to obtain the individual process rates 
is transition-state theory (TST) (Glasston, Laidler and 
Eyring, 1941; Vineyard, 1957; Laidler, 1987). The two 
basic quantities entering this theory are an effective at- 
tempt frequency, To, and the minimum energy barrier 
AE that needs to be overcome for the event to take place, 
i.e. to bring the system from the initial to the final state. 
The atomic configuration corresponding to AE is accord- 
ingly called the transition state (TS). Within a harmonic 
approximation, the effective attempt frequency is pro- 
portional to the ratio of normal vibrational modes at the 
initial and transition state. Just like the barrier AE', 
Fo is thus also related to properties of the PES, and as 
such directly amenable to a calculation with electronic 
structure theory methods like DFT (Ratsch and Schef- 
fler, 1998). 

In the end, the crucial additional PES information re- 
quired in kMC compared to equilibrium MC is therefore 
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FIG. 13: Schematic top view of a fcc(lOO) surface, explaining 
diffusion processes of an isolated metal adatom (white circle), 
a) Diffusion by hopping to a neighboring lattice site, b) dif- 
fusion by exchange with a surface atom. 



the location of the transition state in form of the PES 
saddle point along a reaction path of the process. Par- 
ticularly for high-dimensional PES this is not at all a 
trivial problem, and the development of efficient and re- 
liable transition-state-search algorithms is a very active 
area of current research (Henkelman, Johannesson and 
Jonsson, 2000). For many surface related elementary pro- 
cesses (e.g. diffusion, adsorption, desorption or reaction 
events) the dimensionality is fortunately not excessive, 
or can be mapped onto a couple of prominent reaction 
coordinates as exemplified in Figure 1121 The identifica- 
tion of the TS and the ensuing calculation of the rate for 
individual identified elementary processes with TST are 
then computationally involved, but just feasible. 

This still leaves as a fundamental problem, how the rel- 
evant elementary processes for any given system configu- 
ration can be identified in the first place. Most TS-search 
algorithms require not only the automatically provided 
information of the actual system state, but also knowl- 
edge of the final state after the process has taken place 
(Henkelman, Johannesson and Jonsson, 2000). In other 
words, quite some insight into the physics of the elemen- 
tary process is needed in order to determine its rate and 
include it in the list of possible processes in the kMC 
simulation. How difficult and non-obvious this can be 
even for the simplest kind of processes is nicely exem- 
plified by the diffusion of an isolated metal atom over 
a close-packed surface (Ala-Nissila, Ferrando and Ying, 
2002). Such a process is of fundamental importance for 
the epitaxial growth of metal films, which is a necessary 
prerequisite in many applications like catalysis, magneto- 
optic storage media or interconnects in microelectronics. 
Intuitively, one would expect the surface diffusion to pro- 
ceed by simple hops from one lattice site to a neighboring 
lattice site, as illustrated in Figure 113b . for a fee (100) 
surface. Having said that, it is in the meanwhile well es- 
tablished that on a number of substrates diffusion does 
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not operate preferentially by such hopping processes, but 
by atomic exchange as explained in Figure lCT p. Here, the 
adatom replaces a surface atom, and the latter then as- 
sumes the adsorption site. Even much more complicated, 
correlated exchange diffusion processes involving a larger 
number of surface atoms are currently discussed for some 
materials. And the complexity increases of course fur- 
ther, when diffusion along island edges, across steps and 
around defects needs to be treated in detail (Ala-Nissila, 
Ferrando and Ying, 2002). 

While it is therefore straightforward to say that one 
wants to include e.g. diffusion in a kMC simulation, it 
can in practice be very involved to identify the individual 
processes actually contributing to it. Some attempts 
to automatize the search for the elementary processes 
possible for a given system configuration are currently 
undertaken, but in the large majority of first-principles 
kMC studies performed up to date and in the foreseeable 
future, the process lists are simply generated by physical 
insight. This obviously bears the risk of overlooking a 
potentially relevant molecular process, and on this note 
this just evolving method has to be seen. Contrary to 
traditional kMC studies, where an unknown number 
of real molecular processes is often lumped together 
into a handful effective processes with optimized rates, 
first-principles kMC has the advantage, however, that 
the omission of a relevant elementary process will defi- 
nitely show up in the simulation results. As such, first 
experience tells that a much larger number of molecular 
processes needs to be accounted for in a corresponding 
modeling "with microscopic understanding" compared 
to traditional empirical kMC (Stampfl et ai, 2002). In 
other words, that the statistical interplay determining 
the observable function of materials takes places between 
quite a number of different elementary processes, and 
is therefore often way too complex to be understood 
by just studying in detail the one or other elementary 
process alone. 

Applications to semiconductor growth and catalysis 

The new quality of and the novel insights that can 
be gained by mesoscopic first-principles kMC simula- 
tions was first demonstrated in the area of nucleation 
and growth in metal and semiconductor epitaxy (Rug- 
gerone, Ratsch and Scheffler, 1997; Ratsch, Ruggerone 
and Scheffler, 1998; Ovesson, Bogicevic and Lundqvist, 
1999; Fichthorn and Scheffler, 2000; Kratzer and Schef- 
fler, 2001; Kratzer and Scheffler, 2002; Kratzer, Penev 
and Scheffler, 2002). As one example from this field we 
return to the GaAs(OOl) surface already discussed in the 
context of the chemical potential plots. As apparent from 
Figure^ the so-called /32(2 x 4) reconstruction repre- 
sents the most stable phase under moderately As-rich 
conditions, which are typically employed in the molecu- 
lar beam epitaxy (MBE) growth of this material. Aiming 
at an atomic-scale understanding of this technologically 
most relevant process, first-principles LGH -I- kMC sim- 




FIG. 14: Snapshots of characteristic stages during a first- 
principles kMC simulation of GaAs homoepitaxy. Ga and As 
substrate atoms appear in green and dark blue, Ga adatoms 
in yellow, and freshly adsorbed As dimers in light blue, (a) 
Ga adatoms preferentially wander around in the trenches, (b) 
Under the growth conditions used here, an As2 molecule ad- 
sorbing on a Ga adatom in the trench initiates island forma- 
tion, (c) Growth proceeds into a new atomic layer via Ga 
adatoms forming Ga dimers. (d) Eventually, a new layer of 
arsenic starts to grow, and the island extends itself towards 
the foreground, while more material attaches along the trench 
(from Kratzer and Scheffler, 2002). 



ulations were performed, including the deposition of As2 
and Ga from the gas-phase, as well as diffusion on this 
complex /32(2 x 4) semiconductor surface. In order to 
reach a trustworthy modeling, the consideration of more 
than 30 different elementary processes was found to be 
necessary, underlining our general message that complex 
materials properties cannot be understood by analyzing 
isolated molecular processes alone. Snapshots of char- 
acteristic stages during a typical simulation at realistic 
deposition ffuxes and temperature are given in Figure 
IT^ They show a small part of the total mesoscopic sim- 
ulation area, focusing on one "trench" of the f32{2 x 4) 
reconstruction. At the chosen conditions, island nucle- 
ation is observed in these reconstructed surface trenches, 
which is followed by growth along the trench, thereby 
extending into a new layer. 

Monitoring the density of the nucleated islands in huge 
simulation cells (160 x 320 surface lattice constants), 
a saturation indicating the beginning of steady-state 
growth is only reached after simulation times of the order 
of seconds for quite a range of temperatures. Obviously, 
neither such system sizes, nor time scales would have 
been accessible by direct electronic structure theory cal- 
culations combined e.g. with MD simulations. In the 
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FIG. 15: Saturation island density corresponding to steady- 
state MBE of GaAs as a function of the inverse growth tem- 
perature. The dashed line shows the prediction of classical 
nucleation theory (CNT) for diffusion-limited attachment and 
a critical nucleus size equal to 1. The significant deviation at 
higher temperatures is caused by arsenic losses due to desorp- 
tion, which is not considered in CNT (from Kratzer, Penev 
and Schefller, 2002). 
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FIG. 16: Time evolution of the site occupation by O and 
CO of the two prominent adsorption sites of the RuO2(110) 
model catalyst surface shown in Figure |K| The temperature 
and pressure conditions chosen (T = 600 K, pco ~ 20 atm, 
PO2 = 1 atm) correspond to an optimum catalytic perfor- 
mance. Under these conditions kinetics builds up a steady- 
state surface population in which O and CO compete for ei- 
ther site type at the surface, as refiected by the strong fiuctu- 
ations in the site occupations. Note the extended time scale, 
also for the "induction period" until the steady-state popula- 
tions are reached when starting from a purely oxygen covered 
surface (from Renter, Frenkel and Scheffier, 2004). 



ensuing steady-state growth, attachment of a deposited 
Ga atom to an existing island typically takes place before 
the adatom could take part in a new nucleation event. 
This leads to a very small nucleation rate that is coun- 
terbalanced by a simultaneous decrease in the number of 
islands due to coalescence. The resulting constant island 
density during steady-state growth is plotted in Figure lT^ 
for a range of technologically relevant temperatures. At 
the lower end around 500-600 K, this density decreases, as 
is consistent with the frequently employed classical nucle- 
ation theory (CNT). Under these conditions, the island 
morphology is predominantly determined by Ga surface 
diffusion alone, i.e. it may be understood on the basis 
of one molecular process class. Around 600 K the island 
density becomes almost constant, however, and even in- 
creases again above around 800 K. The determined mag- 
nitude is then orders of magnitude away from the predic- 
tion of CNT, cf. Figure El but in very good agreement 
with existing experimental data. The reason for this un- 
usual behavior is that the adsorption of As2 molecules 
at reactive surface sites becomes reversible at these ele- 
vated temperatures. The initially formed Ga-As-As-Ga2 
complexes required for nucleation, cf. Figure 114b . be- 
come unstable against As2 desorption, and a decreasing 
fraction of them can stabilize into larger aggregates. Due 
to the contribution of the decaying complexes, an effec- 
tively higher density of mobile Ga adatoms results at 
the surface, which in turn yields a higher nucleation rate 
of new islands. The temperature window around 700- 
800 K, which is frequently used by MBE crystal growers, 
may therefore be understood as permitting a compro- 



mise between high Ga adatom mobility and stability of 
As complexes that leads to a low island density and cor- 
respondingly smooth films. 

Exactly under the technologically most relevant con- 
ditions, the desired functionality of the surface results 
therefore from the concerted interdependence of distinct 
molecular processes, i.e. in this case diffusion, adsorption 
and desorption. To further show that this is to our opin- 
ion more the rule than an exception in materials science 
applications, we return in the remainder of this section to 
the field of heterogeneous catalysis. Here, the conversion 
of reactants into products by means of surface chemical 
reactions {A + B ^ C) adds another qualitatively differ- 
ent class of processes to the statistical interplay. In the 
context of the thermodynamic chemical potential plots 
we had already discussed that these on-going catalytic 
reactions at the surface continuously consume the ad- 
sorbed reactants, driving the surface populations away 
from their equilibrium value. If this has a significant ef- 
fect, presumably e.g. in regions of very high catalytic 
activity, the average surface coverage and structure does 
even under steady-state operation never reach its equilib- 
rium with the surrounding reactant gas-phase, and must 
thence be modeled by explicitly accounting for the sur- 
face kinetics (Hansen and Neurock, 1999; Hansen and 
Neurock, 2000; Renter, Frenkel and Scheffler, 2004). 

In terms of kMC, this means that in addition to the 
diffusion, adsorption and desorption of the reactants and 
products, also reaction events have to be considered. For 
the case of CO oxidation as one of the central reactions 
taking place in our car catalytic converters, this trans- 
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FIG. 17: Left panel: Steady state surface structures of 
RuO2(110) in an O2/CO environment obtained by first- 
principles kMC calculations at T = 600 K. In all non- white 
areas, the average site occupation is dominated (> 90%) by 
one species, and the site nomenclature is the same as in Fig- 
ure |3 where the same surface structure was addressed within 
the ab initio atomistic thermodynamics approach. Right 
panel: Map of the corresponding catalytic CO oxidation ac- 
tivity measured as so-called turn-over frequencies (TOFs), i.e. 
CO2 conversion per cm^ and second: White areas have a 
TOF < 10^^cm~^s~^, and each increasing gray level repre- 
sents one order of magnitude higher activity. The highest 
catalytic activity (black region, TOF > 10^^cm~^s~^) is nar- 
rowly concentrated around the phase coexistence region that 
was already suggested by the thermodynamic treatment (from 
Renter, Frenkel and Scheffler, 2004). 



lates into the conversion of adsorbed O and CO into 
CO2. Even for the afore discussed, moderately complex 
model catalyst RuO2(110), again close to 30 elementary 
processes result, comprising both adsorption to and des- 
orption from the two prominent site-types at the surface 
(br and cus, cf. Figure El, as well as diffusion between 
any nearest neighbor site-combination (br— >br, br^cus, 
cus— s-br, cus— >cus). Finally, reaction events account for 
the catalytic activity and are possible whenever O and 
CO are similarly adsorbed in any nearest neighbor site- 
combination. For given temperature and reactant pres- 
sures, the corresponding kMC simulations are then first 
run until steady-state conditions are reached, and the 
average surface populations are thereafter evaluated over 
sufficiently long times. We note that even for elevated 
temperatures, both time periods may again largely ex- 
ceed the time span accessible by current MD techniques 
as exemplified in Figure ITHI The obtained steady-state 
average surface populations at T = 600 K are shown in 
FigurelTTIas a function of the gas-phase partial pressures. 
Comparing with the surface phase diagram of Figure [S] 
from ab initio atomistic thermodynamics, i.e. neglect- 
ing the effect of the on-going catalytic reactions at the 
surface, similarities, but also the expected significant dif- 
ferences under some environmental conditions can be dis- 
cerned. 

The differences affect most prominently the presence 
of oxygen at the br sites, where it is much more strongly 
bound than CO. For the thermodynamic approach only 



the ratio of adsorption to desorption matters, and due to 
the ensuing very low desorption rate, O'^'' is correspond- 
ingly stabilized even when there is much more CO in the 
gas-phase than O2 (left upper part of FigurelS)). The sur- 
face reactions, on the other hand, provide a very efficient 
means of removing this O*^' species that is not accounted 
for in the thermodynamic treatment. As net result, un- 
der most CO-rich conditions in the gas-phase, oxygen is 
faster consumed by the reaction than it can be replen- 
ished from the gas-phase. The kMC simulations covering 
this effect yield then a much lower surface concentration 
of O'^'^, and in turn show a much larger stability range 
of surface structures with CO'^'' at the surface (blue and 
hatched blue regions). It is particularly interesting to no- 
tice, that this yields a stability region of a surface struc- 
ture consisting of only adsorbed CO at br sites that does 
not exist in the thermodynamic phase diagram at all, cf. 
Figure The corresponding CO*^'/- "phase" (hatched 
blue region) is thus a stable structure with defined aver- 
age surface population that is entirely stabilized by the 
kinetics of this open catalytic system. 

These differences were conceptually anticipated in 
the thermodynamic phase diagram, and qualitatively 
delineated by the hatched regions in Figure |S| Due 
to the vicinity to a phase transition and the ensuing 
enhanced dynamics at the surface, these regions were 
also considered as potential candidates for highly ef- 
ficient catalytic activity. This is in fact confirmed by 
the first-principles kMC simulations as shown in the 
right panel of Figure El Since the detailed statistics 
of all elementary processes is explicitly accounted for in 
the latter type simulations, it is straightforward to also 
evaluate the average occurrence of the reaction events 
over long time periods as a measure of the catalytic 
activity. The obtained so-called turnover frequencies 
(TOF, in units of formed CO2 per cm^ per second) are 
indeed narrowly peaked around the phase coexistence 
line, where the kinetics builds up a surface population 
in which O and CO compete for either site type at the 
surface. This competition is in fact nicely reflected by 
the large fluctuations in the surface populations apparent 
in Figure 1161 The partial pressures and temperatures 
corresponding to this high activity "phase", and even 
the absolute TOF values under these conditions, agree 
extremely well with detailed experimental studies mea- 
suring the steady-state activity in the temperature range 
from 300-600 K and both at high pressures and in UHV. 
Interestingly, under the conditions of highest catalytic 
performance it is not the reaction with the highest rate 
(lowest barrier) that dominates the activity. Although 
the particular elementary process itself exhibits very 
suitable properties for catalysis, it occurs too rarely 
in the full concert of all possible events to decisively 
affect the observable macroscopic functionality. This 
emphasizes again the importance of the statistical 
interplay and the novel level of understanding that can 
only be provided by first-principles based mesoscopic 
studies. 
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Outlook 

As highlighted by the few examples from surface 
physics discussed above, many materials' properties and 
functions arise out of the interplay of a large number 
of distinct molecular processes. Theoretical approaches 
aiming at an atomic-scale understanding and predictive 
modeling of such phenomena have therefore to achieve 
both an accurate description of the individual elementary 
processes at the electronic regime and a proper treat- 
ment of how they act together on the mesoscopic level. 
We have sketched the current status and future direction 
of some emerging methods which correspondingly try to 
combine electronic structure theory with concepts from 
statistical mechanics and thermodynamics. The results 
already achieved with these techniques give a clear indi- 
cation of the new quality and novelty of insights that can 
be gained by such descriptions. On the other hand, it is 
also apparent that we are only at the beginning of a suc- 
cessful bridging of the micro- to mesoscopic transition in 
the multi-scale materials modeling endeavor. Some of the 
major conceptual challenges we see at present that need 
to be tackled when applying these schemes to more com- 
plex systems have been touched in this chapter. They 
may be summarized under the keywords accuracy, map- 
ping and efficiency, and as outlook we briefly comment 
further on them. 

Accuracy: The reliability of the statistical treatment 
depends predominantly on the accuracy of the descrip- 
tion of the individual molecular processes that are input 
to it. For the mesoscopic methods themselves it makes 
in fact no difference, whether the underlying PES comes 
from a semi-empirical potential or from first-principles 
calculations, but the predictive power of the obtained re- 
sults (and the physical meaning of the parameters) will 
obviously be significantly different. In this respect we 
only mention two somehow diverging aspects. For the 
interplay of several (possibly competing) molecular pro- 
cesses, an accurate absolute description of each individual 
process e.g. in form of a rate for kMC simulations may 
be less important than the relative ordering among the 
processes as e.g. provided by the correct trend in their 
energetics. In this case, the frequently requested chemi- 
cal accuracy in the description of single processes could 
be a misleading concept, and modest errors in the PES 
would tend to cancel (or compensate each other) in the 
statistical mechanics part. For the particular case of DFT 
as the current workhorse of electronic structure theories 
this could mean that the present uncertainties due to the 
approximate treatment of electronic exchange and corre- 
lation are less problematic than hitherto often assumed. 
On the other hand, in other cases where for example one 
process strongly dominates the concerted interplay, such 
a cancellation will certainly not occur. Then, a more ac- 
curate description of this process will be required than 
can be provided by the exchange-correlation functionals 
in DFT that are available today. Improved descriptions 



based on wave-function methods and on local corrections 
to DFT exist or are being developed, but come so far at 
a high computational cost. Assessing what kind of ac- 
curacy is required for which process under which system 
state, possibly achieved by evolutionary schemes based 
on gradually improving PES descriptions, will therefore 
play a central role in making atomistic statistical me- 
chanics methods computationally feasible for increasingly 
complex systems. 

Mapping: The configuration space of most materials 
science problems is exceedingly large. In order to arrive 
at meaningful statistics, even the most efficient sampling 
of such spaces still requires (at present and in the fore- 
seeable future) a number of PES evaluations that is pro- 
hibitively large to be directly provided by first-principles 
calculations. This problem is mostly circumvented by 
mapping the actual system onto a coarse-grained lat- 
tice model, in which the real Hamiltonian is approxi- 
mated by discretized expansions e.g. in certain interac- 
tions (LGH) or elementary processes (kMC). The expan- 
sions are then first parametrized by the first-principles 
calculations, while the statistical mechanics problem is 
thereafter solved exploiting the fast evaluations of the 
model Hamiltonians. Since in practice these expansions 
can only comprise a finite number of terms, the mapping 
procedure intrinsically bears the problem of overlooking 
a relevant interaction or process. Such an omission can 
obviously jeopardize the validity of the complete statis- 
tical simulation, and there are at present no fool-proof 
or practical, let alone automatized schemes as to which 
terms to include in the expansion, neither how to judge 
on the convergence of the latter. In particular when go- 
ing to more complex systems the present "hand-made" 
expansions that are mostly based on educated guesses 
will become increasingly cumbersome. Eventually, the 
complexity of the system may become so large, that even 
the mapping onto a discretized lattice itself will be prob- 
lematic. Overcoming these limitations may be achieved 
by adaptive, self-refining approaches, and will certainly 
be of paramount importance to ensure the general appli- 
cability of the atomistic statistical techniques. 

Efficiency: Even if an accurate mapping onto a model 
Hamiltonian is achieved, the sampling of the huge con- 
figuration spaces will still put increasing demands on the 
statistical mechanics treatment. In the examples dis- 
cussed above, the actual evaluation of the system par- 
tition function e.g. by (k)MC simulations is a small add- 
on compared to the computational cost of the underlying 
DFT calculations. With increasing system complexity, 
different problems and an increasing number of processes 
this may change eventually, requiring the use of more ef- 
ficient sampling schemes. A major challenge for increas- 
ing efficiency is for example the treatment of processes 
operating at largely different time scales. The computa- 
tional cost of a certain time span in kMC simulations is 
dictated by the fastest process in the system, while the 
slowest process governs what total time period needs ac- 
tually to be covered. If both process scales differ largely. 
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kMC becomes expensive. Remedy may e.g. be provided 
by assuming the fast process to be always equilibrated 
at the time scale of the slow one, and correspondingly an 
appropriate mixing of equilibrium MC with kMC simula- 
tions may significantly increase the efficiency (as already 
done in nowadays TPD simulations). Alternatively, the 
fast process could not be explicitly considered anymore 
on the atomistic level, and only its effect incorporated 
into the remaining processes. 

Obviously, with such a grouping of processes one ap- 
proaches already the meso- to macroscopic transition, 
gradually giving up the atomistic description in favor of 
a more coarse-grained or even continuum modeling. The 
crucial point to note here is that such a transition is done 
in a controlled and hierarchical manner, i.e. necessarily 
as the outcome and understanding from the analysis of 
the statistical interplay at the mesoscopic level. This is 
therefore in marked contrast to e.g. the frequently em- 
ployed rate equation approach in heterogeneous catalysis 
modeling, where macroscopic differential equations are 
directly fed with microscopic parameters. If the latter 
are simply fitted to reproduce some experimental data, 
at best a qualitative description can be achieved anyway. 
If really microscopically meaningful parameters are to be 
used, one does not know which of the many in princi- 
ple possible elementary processes to consider. Simple- 
minded "intuitive" approaches like e.g. parametrizing 
the reaction equation with the data from the reaction 
process with the highest rate may be questionable in 
view of the results described above. This process may 
never occur in the full concert of the other processes, 
or it may only contribute under particular environmen- 
tal conditions, or be significantly enhanced or suppressed 
due to an intricate interplay with another process. All 
this can only be filtered out by the statistical mechanics 
at the mesoscopic level, and can therefore not be grasped 
by the traditional rate equation approach omitting this 
intermediate time and length scale regime. 

The two key features of the atomistic statistical 
schemes reviewed here are in summary that they treat 
the statistical interplay of the possible molecular pro- 
cesses, and that these processes have a well-defined mi- 
croscopic meaning, i.e. they arc described by parameters 
that are provided by first-principles calculations. This 
distinguishes these techniques from approaches where 
molecular process parameters arc cither directly put into 
macroscopic equations neglecting the interplay, or where 
only effective processes with optimized parameters are 
employed in the statistical simulations. In the latter 
case, the individual processes loose their well-defined mi- 
croscopic meaning and typically represent an unspecified 
lump sum of not further resolved processes. Both the 
clear cut microscopic meaning of the individual processes 
and their interplay are, however, decisive for the trans- 
ferability and predictive nature of the obtained results. 
Furthermore, it is also precisely these two ingredients 
that ensure the possibility of reverse-mapping, i.e. the 
unambiguous tracing back of the microscopic origin of 



(appealing) materials' properties identified at the meso- 
or macroscopic modeling level. We are convinced that 
primarily the latter point will be crucial when trying to 
overcome the present trial and error based system en- 
gineering in materials sciences in the near future. An 
advancement based on understanding requires theories 
that straddle various traditional disciplines. The ap- 
proaches discussed here employ methods from various 
areas of electronic structure theory (physics as well as 
chemistry), statistical mechanics, mathematics, materi- 
als science, and computer science. This high interdisci- 
plinarity makes the field challenging, but is also part of 
the reason why it is exciting, timely and full with future 
perspectives. 
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